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ABSTRACT 

Aims. We test cosmological models of structure formation using the rotation curve of the nearest spiral galaxy, M31, determined using 
a recent deep, full-disk 21 -cm imaging survey smoothed to 466 pc resolution. 

Methods. We fit a tilted ring model to the HI data from 8 to 37 kpc and establish conclusively the presence of a dark halo and its 
density distribution via dynamical analysis of the rotation curve. 

Results. The disk of M31 warps from 25 kpc outwards and becomes more inclined with respect to our line of sight. Newtonian 
dynamics without a dark matter halo provide a very poor fit to the rotation curve. In the framework of modified Newtonian dynamic 
(MOND) however the 21-cm rotation curve is well fitted by the gravitational potential traced by the baryonic matter density alone. 
The inclusion of a dark matter halo with a density profile as predicted by hierarchical clustering and structure formation in a ACDM 
cosmology makes the mass model in newtonian dynamic compatible with the rotation curve data. The dark halo concentration param- 
eter for the best fit is C = 12 and its total mass is 1.2 10'" Mq. If a dark halo model with a constant-density core is considered, the 
core radius has to be larger than 20 kpc in order for the model to provide a good fit to the data. We extrapolate the best-fit ACDM and 
constant-density core mass models to very large galactocentric radii, comparable to the size of the dark matter halo. A comparison of 
the predicted mass with the M31 mass determined at such large radii using other dynamical tracers, confirms the validity of our results. 
In particular the ACDM dark halo model which best fits the 21-cm data well reproduces the mass of M31 traced out to 560 kpc. Our 
best estimate for the total mass of M31 is 1.3x10'^ M©, with 12% baryonic fraction and only 6% of the baryons in the neutral gas 
phase. 

Key words. Galaxies: individual (M 31) - Galaxies: ISM - Galaxies: kinematics and dynamics - radio lines: galaxies - dark matter 



^ 1. Introduction ample iNavarro et al.l (|2004') proposed the "Einasto profile", a 

three-parameter formulation, to improve the accuracy of the fits 
j^j Rotation curves of spiral galaxies are fundamental tools to study to cuspy inner density profiles of simulated halos. The two pa- 
d the visible mass distributions in galaxies and to infer the prop- rameters of the "universal" NEW density profile are the halo 
erties of any associated dark matter halos. These can then be overdensity and the scale radius, or (in a more useful param- 
used to constrain cosmological models of galaxy formation and eterization) the halo concentration and its virial mass. Eor hi- 
evolution. Great effort has been devoted in recent years to test erarchical structure formation, small galaxies should show the 
theoretical predictions of cosmological models regarding the de- highest halo concen trations and massive galaxies the lowest ones 
tailed structure of dark matter halos via observations on galac- (iMaccio et al.ll2007 ). Do the observations confirm these predic- 
tic and sub-galactic scales. Knowledge of the halo density pro- tions? Dwarf galaxies with extended rotation curves have often 
file from the center to the outskirts of galaxies is essential for contradicted this theory since central regions show shallo w den- 
solving crucial issues at the heart of galaxy formation theories, sity cores, i.e. very low dark m atter concentrations (.Rhee et al.l 
including the nature of the dark matter itself. Numerical simu- 5004t lGentile"etaLll200l l2007h . This has given new insights on 
lations of structure formation in the flat Cold Dark Matter cos- the nature of dark matter and lead to discussion on how the halo 
mological scenario (hereafter ACDM) predict a well defined ra- structure might have been altere d by the galaxy formation pro- 
dial density profile for the c oUisionless particles in virialized cess (e.g. lOnedin & Zhaoll2002h . Recent hydrodynamical sim- 
structures, the NEW profile (iNavarro et al.l 1 19961) . While there ulations in ACDM framework have shown that strong outflows 
is a general consensus that hierarchical assembly of A CDM in dwarf galaxies inhibit the formation of bulges and decrease 
halos yields "universal" mass profiles on a large scale (i.e. in- the dark matter density, thus reconciling dwarf galaxies with 
dependent of mass and cosmology aside from a simple two pa- ACDM theoretical predictions (Governato et al. 2009). Shallow 
rameter scaling), there is still some controversy on the central density cores have often been supported also by high resolution 
density profile, and on the relative scaling parameters. Eor ex- 
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in today's galaxies dvan den B osch et al. 2000; Corbelli 2003; 



analysis of rotation curves of spirals and low sur face brightness 
galaxies dde Blok et al.ll2001l:lGentile et al.ll2004 '). However, un- 
certainties related to the presence of non-circular motion (often 
related to the presence of small bars and bulges), observational 
uncertainties on the observed velocities, and the possibility of 
dark matter compression during the baryonic collapse leave still 
open the question o n how d ark matter is effectively distributed 

( van den 

ICorbelU & Walterbosll2007ir 

Bright galaxies, because of the large fraction of visible to 
dark matter, do not offer the possibility to trace dark matter very 
accurately in the center. Uncertainties related to distance esti- 
mates, to the disk-bulge light decompositions, and the typically 
limited extention of the gaseous disks beyond the bright star 
forming disks further limit the ability to derive accurate dynam- 
ical mass models. These difficulties can be alleviated in the case 
of the spiral galaxy M31 (Andromeda). Owing to its size, prox- 
imity, well known distance, and to constraints on its structural 
parameters from the long history of observations at all wave- 
lengths, Andromeda (the nearest giant spiral galaxy) offers a 
unique opportunity to analyze in detail the mass distribution and 
the dark halo properties in bright disk galaxies. Massive galaxies 
like M3 1 can probe dark matter on mass scales much larger than 
that of the dwaifs, of order 10'^ Mq. 

The Milky Way and Andromeda are the most massive mem- 
bers of the Local Group. Any estimate of their total mass and of 
the structure of their dark matter halos is a requirement for any 
study of the dynamics of the Local Group, i ts forma tion, evolu- 
tion, and ultimate fate of its members (e.g. iLi & W hite 2008). 
Difficulties in the determination of the Milky Way's mass com- 
ponents, related to the fact that our solar system is deeply em- 
bedded in its disk, can be overcome in the case of M31. M31 is 
known to have a complex merging history. Its multiple nucleus 
(Lauer et al. 1993 ) and the extended stel lar stream and halo (e.g. 
flrwin et al.1120011: IChapman et al.ll2008l) are clear signs of a tu- 
multuous life. According to the hierarchical models of galaxy 
formation it is conceivable that M3 1 has grown by accretion of 
numerous small galaxies. It is likely the most massive member 
of the Local Group (e.g. Klypin et al. 2002). It is therefore of 
great interest to test the other predictions of hierarchical models 
such as the presence and structure of a dark matter halo around 
it. Contrary to dwarf galaxies, luminous high surface brightness 
galaxies such as M31, cannot be used to test the central dark 
matter distribution, not only because of the large surface den- 
sity of baryons which makes it difficult to constrain dark matter, 
but also because o f possible adiabatic contraction ( Klvpin et al.l 
120021: ISeigar et al. 112008). An extended and well defined rotation 
curve can instead be complemented by the extensive informa- 
tion now available on the M3 1 stellar disk, stellar stream, globu- 
lar clusters and orbits of Andromeda's small satellite galaxies to 
establish the dark matter density profile at large galactocentric 
distances. And this is one of our goals. 

M31 was one of the first galaxies where Slipher (in 1914) 
found evidence of rotati on and also the first g alaxy to have a 
published velocity field dSofue & RubinI 200 iL and re ferences 
therein). Using the M31 rotation curve. iBabcockl (1 1 939h was the 
first person to advocate unseen mass at large radii in a galaxy. 
Since then much effort has been devoted to study the rotation 
curve of M3 1 and to understand the relation between the light 
and the mass distribution. Despite a century of dedicated work, 
there are still many unsettled questions concerning the shape of 
the M3 1 rotation curve, the contribution of visible and dark mat- 
ter to it, and the changing orientation of the M3 1 disk. Detailed 
HI surveys with single dish or synthesis observations (e.g. 



Newton & ErnersonI 119771: ICram et al.1 119801: iBaiaia & Shaiig 
il982t lUnwinlll983t [Brinks & Burtonlll984t lBraunlll991l) have 



been analyzed to find local kinematical signatures of spiral arm 
segments, of the interaction with M32 or of a warped disk or 
ring. Even though many authors have pointed out the presence 
of a warp in M31, i.e. of a systematic deviation of the matter 
distribution from equatorial symmetry, a complete quantitative 
analysis of the parameters of such a distortion at large galacto- 
centric radii is still missing. Previous modelling of the HI warp 
has been based on a combination of high resolution inner disk 
HI data and much lower resolution and sensitivity outer disk 
data. The models assumed a rotation curve but no independent 
fit of rotation curv e and warp ing of the disk has been attempted 



(e.g. Henderso"nlll 979; Bri nks & BurtonI 11984 . Some papers 



(e.g. ICarignan et al.ii2006) analyze the extended rotation curve 
of M3 1 using only HI data along the direction of the optical ma- 
jor axis, without c onsidering the poss ibility of a warped disk. 
Only very recently IChemin et al.l (12009) use deep 21 -cm survey 
of the M31 based on high resolution synthesis observations to 
model the warp and the rotation curve simultaneously. 

Our first a im is to use the new WRST HI survey of M31 
(iBraun et al.l2009.) to define the amplitude and orientation of the 
warp using a tilted ring model. A set of free rings will be consid- 
ered for which the following parameters need to be determined: 
the orbital center, the systemic velocity, the inclination and po- 
sition angle with respect to our line of sight, and the rotational 
velocity. The geometric properties of the best fitting tilted ring 
model will then be used to derive the rotation curve from the 
21 -cm line observed velocities. The final goal will be to use the 
rotation curve for constraining the baryonic content of the M3 1 
disk and the presence and distribution of dark matter in its halo 
through the dynamical analysis. 

Our recent deep wide-field HI imaging survey reaches a 
maximum resolution of about 50 pc and 2 km s"' across a 
95x48 kpc^ region. This makes our database the most detailed 
ever made of the neutral medium of any complete galaxy disk, 
including our own. Observations and data reduction are de- 
scribed in IBraun et al.l (1200^) (hereafter Paper I). In Paper I we 
analyzed HI self-absorption features and find opaque atomic gas 
organized into filamentary complexes. While the gas is not the 
dominant baryonic component in M31, we take these opacity 
corrections into account in determining the dynamical contribu- 
tions of the various mass components to the M3 1 rotation curve. 
In this paper we use the data at a resolution of 2 arcmin (457 kpc) 
in order to gain sensitivity in the outermost regions. At this 
spatial resolution, we reach a brightness sensitivity of 0.25 K. 
Considering a typical signal width of 20 km s ' our sensitivity 
should be appropriate for detecting HI gas at column densities 
as low as 10''' cm"^. 

In Section 2, we describe the modeling procedures for deter- 
mining the disk warp in M3 1 and discuss the resulting disk ori- 
entation. In Section 3, we determine the rotation curve and the 
uncertainties associated with it. Various dynamical mass models 
for the rotation curve fit are introduced and discussed in Section 
4. We determine the total baryonic and dark mass of this galaxy. 
Together with complementary data at very large galactocentric 
radii, we confirm the predictions of ACDM cosmological mod- 
els. Section 5 summarizes the main results of this paper. 

We a ssume a distance to M31 o f 785 kpc throughout, as de- 
rived bv iMcConnachie et al.l (l2005h (D=785 + 25 kpc). 
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Fig. 1. The first moment map. The intensity-weighted mean ve- 
locity has been computed from the 120 arcsec data cube at a 
spectral resolution of 2 km s"', using a 4-cr blanking in the data 
cube at 18.5 km/s resolution to determine what is included in the 
integral. 



2. Tilted rings: modeling procedures 

For a dynamical mass model of a disk galaxy it is necessary 
to reconstruct the tri-dimensional velocity field from the veloc- 
ities observed along the Une of sight. If velocities are circular 
and confined to a disk one needs to establish the disk orienta- 
tion for deriving the rotation curve, i.e. the position angle of the 
major axis (PA), and the inclination of the disk with respect to 
the line of sight (/). If the disk exhibits a warp these parame- 
ters vary with galactocentric radius. This is often the case for 
gaseous disks which extend outside the optical radius and which 
often show a different orientation than the inner one. Our at- 
tempt to understand the kinematics of M31 is done performing 
a tilted ring model fit to the data, under the assumption of cir- 
cular motion. Because of this assumption we will use the tilted 
ring model outside the inner 8 kpc region, i.e. where deviations 
of gas motion from circular orbits are expected to be small. We 
will not consider local perturbations to the circular velocity field 
such as those due to spiral arms. The comparison between the 
velocities predicted by a tilted ring model and the data is done 
over all azimuthal angles and not only around the major axis. 
This will average out spiral arm perturbations. 

In Figure [1] we show the first moment map (i.e. the intensity- 
weighted mean velocity along the line of s ight) of our 21 -cm 
Andromeda survey (see iBraun et al.ll2009l for more details). 
Because of the large angular extent of Andromeda, it is neces- 
sary to consider the correct transformation between angular and 
cartesian coordinates in order to derive galactocentric distances. 
A detailed description of the spherical trigonometry involved can 
be found in the literature (e.g. van der Marel & Cioni 2001). We 
shall summarize here the most important formulae. Note that 
these take into account the variations in the distance between us 
and the Andromeda regions that do not lie along the line of nodes 
(the disk on the near side (West) half is closer to us than on the 
far side (East) half). Consider a point at a given right ascension 
and declination (a, S) in the Andromeda disk having inclination 
i and position angle P.A. - 0. The distance between the observer 



and this point is D. The center of the galaxy has right ascension, 
declination (ao, 6q) and distance Do from the observer. Consider 
<p as the angle between the tangent to the great circle on the celes- 
tial sphere through (a, 6) and (ao, <5o) and the circle of constant 
declination 6o measured counterclockwise starting from the axis 
that runs in the direction of decreasing R.A. at constant declina- 
tion 6o (in practice (p = P.A. + 90°). The value of (p along the line 
of nodes (for points along the major axis) is (© = -i- 90°). 
The angular distance between (a, 6) and (ao, ^o) is the second 
angular coordinate called p. We shall work in a Cartesian coor- 
dinate system (x, y, z) which has its origin at the galaxy center, 
the X-axis antiparallel to R.A., the y-axis parallel to declination 
axis and the z-axis toward the observer. The coordinates of the 
observed pixels in this system will be: 

X - D sinp cos(f> y = D sinp sm<p (1) 
where 

D - Dq cos/ /[cos/ cosp - sin/ sinp sin(0 - 0)] (2) 

Because of the warp, / and 9 are not constant. This implies that 
for a given ring model the numerical code starts with a guessed 
matrix for / and a guessed matrix for 6 over the observed pixels. 
Once galactocentric distances are determined it checks the ma- 
trices with the tilted ring model and iterates. Positions and ve- 
locities of a ring at a given inclination and position angle (/,, 0,) 
can be easily derived into this reference system by a two angle 
rotation (/,-, 0,) of the Cartesian coordinate system which has the 
ring in the {x' ,y') plane (e.g. Begeman 1989). 

The tilted ring model is (as usual) based on a number of ge- 
ometrical and kinematical parameters relative to the HI disk. To 
represent the HI distribution in M31 we consider 110 rings be- 
tween 10 and 35 kpc in radius. Each ring is characterized by its 
radius R and by 6 additional parameters: the circular velocity Vc, 
the inclination /, the position angle 9, the systemic velocity V,,., 
and the position shifts of the orbital center with respect to the 
galaxy center (Axc, Aje). These last 3 parameters allow the rings 
to be non concentric and the gas at each radius to have a velocity 
component perpendicular to the disk (such as that given by gas 
outflowing or infalling into the disk). If the average value of this 
velocity component is uniform and non zero over the whole disk 
it implies that effectively the systemic velocity determined from 
the gas velocity field is different from the assumed one. We shall 
also consider the case of a tilted ring model with uniform val- 
ues of Vsys, Axc, Aje. To test the model each ring is subdivided 
into segments of equal area smaller than the spatial resolution of 
the dataset in use. We convolve the emission from various pieces 
with the beam pattern at each location to compute the velocity 
along the line of sight and its coordinates in a rest frame, defined 
above. In the same rest frame we consider the 21 -cm dataset. 
The telescope beam shape for the dataset we use (a smoothed 
version of the original data) can be well described by a Gaussian 
function with full width half maximum (FWHM) of 2 arcmin. 
With this method we naturally account for the possibility that the 
intensity-weighted mean velocity in a pixel might result from the 
superposition along the line of sight of emission from different 
rings. For each position in the map we compare the observed ve- 
locity y** with the velocity predicted by the model y™'"'. The 
observed velocity is the mean velocity of the HI gas at each po- 
sition. The assignment of a measure of the goodness of the fit 
in this modeling procedure is done using a test on the differ- 
ence between V"^ ' and The noise is uniform in our map 
(see iBraun et al.l l2009) and we assign uniform uncertainties cr,- 
to the observed values V°^\ In order to keep the model sensitive 
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to variations of parameters of the outermost rings, each pixel in 
the map is assigned equal weight. Pixels with higher or lower 21- 
cm surface brightness contribute equally to determine the global 
goodness of the model fit. Since the velocity channel width of 
our survey is about 2 km s"', we arbitrarily set cr,- equal to the 
width of 3 channels (6 km s"'). This is simply a scaling factor 
The equation below defines the reduced x'' which we will use 
throughout this paper 



1 ^ 



(3) 



In our definition of N is the number of positions with HI col- 
umn density along the line of sight greater than 5 x 10'^ cm"^ 
(A^ - 31391) and v is the number of degrees of freedom. It is 
unreasonable to consider 110 free rings i.e. each ring with 6 
degrees of freedom, not only because of the large amount of 
computer time required to find a solution but also because so- 
lutions with so many degrees of freedom are not very robust 
(parameter variations for 1 out of 110 rings does not give sen- 
sible variations to the x^)- On the other hand the use of func- 
tional forms that describe the 6 free parameters as a function 
of R, each with 2-3 free parameters, often do n ot give sat- 
isfactory results (e.g. Cor belli & Schneideiill997l) . Following 
[Corbelli & Schneider (1997) we constructed a model in which 
the properties of 1 1 equally spaced rings could be varied inde- 
pendently. We shall call these the 'free rings'. The first free ring 
is centered at 10 kpc, the last one at 35 kpc. This is because plac- 
ing the first/last free ring at smaller/larger galactocentric distance 
make them not very stable. We extrapolate the geometrical pa- 
rameters of the first/last free ring for 2 kpc inward/outward when 
using the tilted ring model to derive the rotation curve at radii be- 
tween 8 and 37 kpc. The parameters of the 10 rings between 2 
free rings are found by linear interpolation. Our procedure is to 
search for a minimum^^ value considering the parameters of the 
1 10 rings. 

Given the numerous free parameters (v - 66) it is still un- 
reasonable to search for the minimal solution scanning a multi- 
dimensional grid (of 66 dimensions). We use two procedures to 
determine the best fitting model. In the first procedure (hereafter 
PI) we apply the technique of partial minima. We evaluate by 
varying each parameter separately and interpolating to estimate 
the value for a minimum for each parameter We then repeat 
the procedure over a smaller parameter interval around the new 
solution. We start the minimization with the following initial set 
of free parameters: / = 77°, Q = 38°, V = -300 km s^', 
Axc - Ajc - for all free rings. As an initial guess of we 
give the average rotational velocity derived in each free ring by 
deprojecting the observed velocities within a 20° cone of the op- 
tical line of nodes, B = 38° (using ; = 77° and Q = 38° for the 
deprojection). We carried out several optimization attempts un- 
der a variety of initial conditions to avoid that our partial minima 
technique carries the solution towards a relative minimum. 

The best fitting ring model is shown in Figure|2]by the heavy 
continuous line. It gives Ilx^ - 2.18. We display only 3 of the 6 
free parameters for each free ring, namely the inclination /, the 
position angle Q and the systemic velocity Ys\s- The shifts of the 
orbital centers are small and shown in Figure|5] The minimum;\f^ 
is obtained for an inclination which radially decreases by a few 
degrees between 10 and 25 kpc in radius and then increases out 
to 85 ° for the outermost ring. The position angle shows marginal 
variations out to 25 kpc, then it decreases by about 10 ° in the 
outskirts. However, while from 10 to 28 kpc \(m values are 




Fig. 2. Assortment of "free ring" models fit to the HI data accord- 
ing to modeling procedure PI. We display 33 of the 66 free pa- 
rameters: the inclination /, position angle Q and systemic veloc- 
ity Ksv.! of the 1 1 free rings. The heavier continuous lines show 
the best fitting model which gives the minimum solution (cir- 
cled points). For this case we vary one parameter of a given ring 
keeping all others at the best fitting values. Dashed lines indi- 
cate some other combinations of free parameters which give a 



obtained for combinations of / and 6 not very diff'erent than the 
best fitting model, beyond 28 kpc we can find tilted ring mod- 
els which give acceptable fits (with a x^ value within 20% of 
the minimum) with different combinations of / and Q. The short 
dashed line in Figure |2] shows one model for which the inclina- 
tion of the outermost rings does not increase and their position 
angle instead varies between 30° and 45°. Clearly the outermost 
rings are not very stable, due to a lack of 21 -cm emission at 
these large radii. Hence, some care is required when the model 
is extrapolated to even larger radii. The single constant value of 
systemic velocity which minimizes the;^'^ is -306 km s"'. The 
errorbars displayed for the best fitting model indicate parame- 
ter variations for which the x^ of the model is within 20% of the 
minimum. This is done varying just one parameter at a time with 
respect to the combination of parameters which gives the mini- 
mum The rotational velocities obtained from the fit will 
be discussed in the next Section where we compare them to the 
rotation curve of Andromeda derived from the first moment map 
using the geometrical parameters and V,,., of the best tilted ring 
model. 

In a second procedure (hereafter P2) we searched for a min- 
imum by neglecting the radial variations of Vjij, jCcJc. We set 



Ax^ 



Aye 



and y„ 



-306 km s . We have determined 



this value of Ysys as described below. In M3 1 , the integrated 21- 
cm profile shows that the North-East receding side has more gas 
than the South- West approaching side (their ratio being 1.13), 
at intermediate velocities between the central one and the veloc- 
ity of maximum emission. This is likely due to the disturbance 
caused by interaction with M32. Hence we cannot simply com- 
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Fig. 3. Best fitting P2 solution for / and Q obtained when Vj,,, — 
-306 km s"', Axc = 0, Ay^ - 0. All combinations of the 12 
free parameters for the outermost 4 free rings have been consid- 
ered. For the innermost 7 free rings instead the minimum values 
have been obtained by varying parameters of one ring at a time. 
Errorbars indicate parameter variations for which the value 
of the model is within 20% of the minimum. The lower panel 
shows the values of V,,,, determined by the first order moment 
map (see text for details) between 9 and 25 kpc, whose average 
value is -306 km s"' . 

pute the systemic velocity by averaging the observed mean ve- 
locities and weighting these with the 21 -cm surface brightness. 
Considering the galaxy disk between 10 and 25 kpc with an av- 
erage inclination of 77.7° and position angle 38° we compute the 
average systemic velocity in rings. We weight the data with the 
surface brightness intensity. For points south of the minor axis 
we multiply their weight by the ratio of the HI mass in the north- 
ern side to that of the southern side in that the ring. The resulting 
Vjys is shown in Figure[3] Averaging the systemic velocities over 
aU rings we thus obtain a value of -306 km s We searched for 
a minimal solution scanning a multidimensional grid (12 dimen- 
sions) corresponding to parameter variations of the outermost 4 
free rings. The initial grid of variations for the inclination and 
position angle considered are +15°, +10°, +5° around the stan- 
dard values of the inner disk (/ - 11 .6°, - 38°). The maximum 
velocity shifts considered with respect to the initial values of Vc 
are +25, 15, 10,5 km s"'. We then consider finer grids for the 
outermost 4 rings. We have compared the observed velocities to 
the velocities predicted by tilted ring models using all possible 
combinations of parameters for the outermost 4 free rings. For 
the 8 innermost free rings, parameter variations have been con- 
sider only in a 3-dimensional space, namely the minimum;^'^ has 
been found for one ring at a time, considering all combinations 
of /, 6, Vc for that particular ring. As shown in Figure [3] the re- 
sulting free parameters for the best fitting tilted ring model are 
very similar to those obtained using the first method. 
We define the residual velocities as 



V, 



mod 



V, 



ohs 



(4) 



where Vmod are the circular velocities Vc, as parameterized 
in the tilted ring model, projected along the line of sight and 
convolved with the beam function. The observed velocities Vobs 
are the mean velocities detected along the line of sight. 

Residual velocities are generally smaller than 10 km s ' and 
no characteristic patterns are visible across the galaxy. It is worth 
noticing that, if we trace the major axis position angle by a close 
inspection of the first moment map, the position angle seems 
about constant at ~ 38° out to galactocentric distances of 28 kpc, 
then it decreases to about 32° at 32 kpc and finally it increases 
again back to 38° at 38 kpc. Given the consistency between the 
best fitting model for the first and second procedure (Figure |2] 
and Figure [3]) and also between these tilted ring models and the 
inspection of the database, we will not consider models in which 
the position angle increases much above 38° for the outermost 
rings in the rest of this paper. 

2.1. The NEMO results 

To check our results for the best fitting tilted ring model we also 
used the s tandar d least-square fitting technique developed by 
iBegemaiil (Il987h . as impleme nted in the RO TCUR task within 
the NEMO software package dTeubenll 19951) (hereafter P3). The 
galactic disk is subdivided into rings, each of which is described 
by the usual 6 parameters. Starting with the initial estimate 
for the fitting parameters, these are then adjusted iteratively for 
each ring independently until convergence is achieved. We run 
NEMO considering or neglecting the variations of the orbital 
centers and systemic velocity for each ring. The ROTCUR task 
gives less accurate results than our PI method since it does not 
take into account the distance variations between the observer 
and the far/near side of the galaxy and minimizes the free param- 
eters of one ring at a time without subsequent iterations. We ran 
ROTCUR with 24 free rings from 8 to 34 kpc. For the last ring 
the solutions did not converge. We show in Figure|4]the resulting 
/, 9 and Vsys for 3 cases. In one case the fit is done over all data 
points with uniform weight, while varying Vsvs and the orbital 
centers Xc andy^.. In a second attempt, we fitted the data weighted 
with the cosine function of the galactic angle away from the ma- 
jor axis and excluding data within a 20° angle around minor axis. 
In a third attempt, we fitted the data without considering possible 
orbital center shifts and setting Vsys to the average value found by 
ROTCUR when Vsys was allowed to vary. This method confirms 
that the average systemic velocity of the disk of Andromeda is 
slightly more negative than previously thought and that there is 
a warp in outer disk which brings the orbits closer to be edge-on. 

We conclude that the 3 fitting methods largely produce sim- 
ilar results and will make use of the results of most accurate 
procedure, PI, in the rest of the paper. 



2.2. Comparison with previous resuits 

As we mentioned in the Introduction, only recently a 21- 
cm survey made with a synthesis telescope has been used to 
model the warp and the rotation curve of M3 1 simultaneously 
(Chemin et al. 2009). The conclusion is similar to ours, namely 
that the warp of M3 1 is such that the outer disk is at higher incli- 
nation and lower position angle than the inner disk. However 
a close inspection reveals some differences. First, the aver- 
age value of the disk inclination between 10 and 28 kpc that 
IChemin et al. (2PC)9) find is 74° while what we find is 77°, more 
consistent with that derived from optical surface photometry 
IWalterbos & Kennicutll (Il987h . Secondly, the inclination of the 
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Fig. 4. Best fitting solutions obtained with the task ROTCUR in 
the NEMO package (P3). The resulting position shifts of orbital 
centers are very small and hence are not shown. The open circles 
indicate the best fitting parameters when we exclude data within 
a 20° angle around minor axis in the plane of the galaxy and the 
data is weighted into the least squares solution with the cosine of 
the galactic angle away from the major axis. The triangles indi- 
cates the solution computed for all data (including points along 
minor axis) and for uniform weight. The star symbols show the 
best fitting solution obtained by keeping the orbital centers fixed, 
and the systemic velocity equal to -305 km s"' , the average value 
obtained when we allow it to vary radially. 



outermost fitted radiu s is higher for our model compared to what 
IChemin et al.l (l2009h find, while the position angle is slightly 
lower. We list below three main differences in the two fitting 
methods which can help explaining the variations in the result- 
ing tilted ring model parameters. 

-We model the bulk rotational velocity of M31 using the 
mean rotational velocity observed at each position. That is, we 
use the first moment map made from the 120 arcsec by 2 km s"' 
data cube, following 4-cr blanking in the 120 arcsec by 18.5 
km s~' data cube to determine what is included in the integral. 
IChemin et all (l2009l) use the peak velocity and if there are mul- 
tiple components they use the peak velocity of the 'main com- 
ponent'. The main component is defined to be the one which 
has the largest velocity relative to the systemic velocity of the 
galaxy. Their choice is based on the observational evidence of 
multiple peak profiles, very prominent in the inner regions. As 
we explain fully in the next subsection we did not include the 
innermost region (from to 8 kpc in radius) in our analysis 
because it is dominated by non-circular motions. Even though 
taking the mean velocity we might be systematically biased to 
lower apparent rotational velocities, the bias effects are mini- 
mal at large radii where we are most interested to the dynamics. 
In general, we find that assuming the peak velocity of the main 
co mponent as the best approximation to the rotation curve, as 
in IChemin et al. I (f2009h . is very model dependent. Moreover, to 
determine it in a robust way one has to assume certain observa- 



Fig. 5. Orbital center shifts in arcmin as derived from tilted 
ring model fits. Positive values of Ax^ correspond to orbital 
centers shifted towards decreasing R.A. values with respect to 
10.6846853° which is Andromeda's center. Positive values of 
Aye correspond to orbital centers shifted towards increasing Dec. 
values with respect to the central value of 41.2690374°. In pan- 
els (fl) and (b) we show the values corresponding to tilted ring 
models PI in Figure|2]where we have considered these variables. 
Similarly in panels (c) and (d) the values derived for models P3 
shown in Figure |4] 



tional conditions (i.e. to which level is a high velocity feature 
accepted in terms of signal-to-noise; also, what happens if the 
highest velocity feature is blended inside a bright component so 
that only a 'shoulder' but not a secondary peak can be seen in the 
spectra). Even though it is true that observations carried out with 
a low resolution and projection effects bias the mean velocity to- 
wards lower velocities, we feel that results might be more robust 
when considering the mean velocity than other estimates if the 
system is as complex, asymmetric and disordered as M3 1 (and 
observed with a high spatial resolution). Accretion events, which 
M3 1 is experiencing, produce several distinctive morphological 
and dynamical signatures in the disk, including long-lived ring- 
like features, significant flares, bars and faint filamentary struc- 
tures above the disk plane. In M3 1 , the likely non-negligible disk 
thickness coupled with a complex, asymmetric warping and with 
the presence of non-circular motion related to multiple spiral 
arm segments intersecting along the line of sight, makes difficult 
to assess the reliability of velocity indicators different than the 
mean for tracing the bulk circular motion of the disk. To prove 
the complexity of the system, we would like to point out that 
across the disk of M3 1 and especially along the bright ring-like 
structure at 10-15 kpc in radius, double peak profiles are often 
present. Sometime peaks are separated by more than 100 km s"' . 
There is however not a systematic azimuthal or radial pattern as 
to whether is the fainter or the brighter peak to show the most ex- 
treme velocity from systemic. If it is just the warped outer disk 
superimposed with the main disk to cause the multiple peak pro- 
files we should have found a more systematic behavior since the 
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neutral hydrogen surface density decreases considerably beyond 
15 kpc. High resolution IR maps of M31 (iGordon et al.ll2006l) 
show that the northern half the 15 kpc arm, distinct on the major 
axis, merges into other spiral arms (or ring like structure) at 10 
kpc. Peak velocities might more closely trace non-circular mo- 
tion of the arms and produce wiggles in the rotation curve which 
do not average out when additional perturbations are present. 
This is the case for M31, in which the southern half is more 
strongly tidally perturbed than the northern half. Such a curve 
cannot be reproduced without modelling the spiral arms locally 
in term of mass condensation and non-circular motion. 

-As we mentioned earlier in this Section we believe that it is 
relevant for M31, a very extended and nearby galaxy, to derive 
galactocentric distances in the frame of a tilted ring model using 
appropriate spherical geometry. This takes into account the fact 
that the near side of the galaxy is effectively closer than the far 
side of the galaxy. Procedures often available for deriving the 
kinematical parameters are built for the more numerous more 
distant galaxies and do not account for this effect. Moreover, 
our routine works with 66 free parameters simultaneously since 
the best fitting values of the parameters for each ring cannot be 
found independently from those relative to other rings when a 
warp is present. In fact the warp makes two or more rings overlap 
along the line of sight and in the overlapping region the expected 
velocity depend on all parameters of the rings which overlap. 

-In fitting the tilted ring model to the data we do not apply 
any angular dependent weight i.e. the same weight is given to 
pixels close to minor axis than to major axis. This is because we 
would like to minimize the risk of amplifying the kinematical 
signatures of sporadic features which happens to lie close to ma- 
jor axis. As we will see in the next Section, our model produces 
a stable rotation curve, not very sensitive to the choice of the 
opening angle a,„ax around the major axis. 

2.3. Why exclude the innermost region from fitting 

Opposite to IChemin et al.l (l2009l) we do not include the in- 
ner regio ns in our analys is. This is because after the pioneer 
work of iLindbladI (Il956h several other papers have pointed 
out in the inner region of M3 1 the presence of morphological 
structures, such as a bar and a bulge, associated with stream- 
ing motion and non -circular o rbits (e.g. f Stark & Binnev 1994; 
Athanassoula & Beaton 2006; Gordon et al. 2006; Beaton et al. 
2007h . In particular Be rman (2001.); Berman & Loinard (.2002) 
have shown that the anomalous velocities observed in the inner 
region of M3 1 can be explained as the response of the gas to 
the potential of a triaxial rotating bulge. Using a bulge effec- 
tive radius of 10 arcmin they have derived which family of peri- 
odic elliptical orbits exist. They find that the bulge gives a non- 
negligible contribution to the galaxy potential out to about 7 kpc, 
and only at larger radii circular motion related to the disk gravi- 
tational potential dominates. The model well reproduced the ve- 
locities observed through the CO J=l-0 line emission. Since in 
this paper we are only modelling the large scale circular motion 
of the disk we use the mean velocity of the HI gas as tracer of 
the circular velocity from 8 kpc outwards. 

3. The rotation curve and the radial distribution of 
the baryons 

We now apply the geometrical parameters of the best fitting tilted 
ring model to derive the rotation curve of the galaxy. We set the 
inclination, position angle and systemic velocity to the values 



shown by the heavy continuous line in Figure |2] and consider 
also the small shifts of the orbital centers obtained by our mini- 
mization procedure PI. We derive the rotation curve by averag- 
ing the rotational velocities of data points in radial bins 1 kpc 
wide. For radii between 8 and 10 kpc, we extrapolate the model 
parameters of the innermost ring centered at 10 kpc. For radii 
larger than 35 kpc, we extrapolate the model parameters of the 
outermost ring at 35 kpc. Just for curiosity, we also checked 
what would be the rotation curve of the inner regions for our 
dataset if we assume circular motion and inside 8 kpc we set 
the disk inclination equal to the inclination of the ring at 10 kpc 
and we consider 28° as position angle (this is the average value 
we derive from an inspection of the mom-1 map). For this in- 
nermost region we consider zero shifts for the orbit centers and 
-306 km s"' as systemic velocity. Figure |6] and Figure |7] show 
the large dispersion in the velocities relative the central region 
due to the presence of multiple components and to the uncer- 
tainties related to orbital eccentricities inside 8 kpc, discussed in 
the previous Section. To complement 21 -cm data in the inner re- 
gions we show in Figure |7] the peak brightness velocities of CO 
lines (iLoinard et al.ll 199511 These have been observed along the 
major axis assumed to be at PA=38° and with ; = 77°. However, 
notice that this is shown not just to point out the consistency of 
the molecular and atomic gas velocities, but to emphasize one 
has to consider non-circular motion to properly tra ce the rota- 
t ion c urve in the inner region (Berman 2001; Berman ^ LoinardI 
12003) . Hence we shall not use the CO data as well as the 21- 
cm data at radii smaller than 8 kpc. In the rest of the paper, 
we will only analyze the rotation curve between 8 and 37 kpc. 
Beyond 37 kpc the northern and southern halves do not give con- 
sistent rotation curves for any of the deconvolution models. This 
is likely due to highly perturbed orbits in the outermost regions, 
especially for the southern half which is closer to M32. 

We consider points which lie within an opening angle a,„ax 
on either side of the major axis. We first derive the rotation curve 
in the northern and southern halves separately. We check the 
consistency of the rotation curves for different values of Qmax, 
with the value of the rotational velocities determined by the tilted 
ring model over the whole galaxy, and between the northern and 
the southern halves. When we vary a„,a^ between 15° and 75°, 
we obtain the same rotation curves consistently (variations are 
less than 1 km s '). Only the dispersion around the mean in- 
creases slightly as we increase a^ax- We shall use Qf,„av - 30° 
for the rest of the paper. The rotation velocities in the two halves 
are consistent (within 3-cr). At many radii, mean velocities in 
the northern and southern halves differ by less than 1-cr corre- 
sponding to only a few km s~' at most, as shown in Figures |6] 
and|7] In Figure|6]we also display the rotational velocity parame- 
ter of the best tilted ring model. The best fitting ring model gives 
the most consistent rotation curves. We also tried to deconvolve 
the observed velocities using other tilted ring models shown in 
Figure |2] but they give less consistent rotation curves. Average 
velocities are derived in the two halves by assigning a weight to 
each pixel equal to the integrated brightness intensity i.e. to the 
HI column density along the line of sight in the limit of optically 
thin 21 -cm line. The global rotation curve is the arithmetic mean 
of the average rotational velocities in the two galaxy halves. The 
errorbars of the global rotation curve shown in the upper panel 
of Figure [T] are computed as 



where Vhi refers to the highest rotational velocity and Vioh- to 
the lowest rotational velocity between the two mean velocities 
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Fig. 6. Rotational velocities derived for the northern and south- 
ern halves of M31 using the geometrical parameters of the best 
fitting tilted ring model as shown in Figure |2] For this Figure 
Q'mflA = 30°. The filled circles indicate the best fit rotational ve- 
locity parameters derived by the tilted ring model fit PI . 



relative to the northern and to the southern galaxy side; cr/„ and 
(Tiotv refer to the relative dispersions around the mean. 

In Table 1 we give the parameters of the best fitting tilted 
ring model and of the average rotation curve of M3 1 in the radial 
interval which will be used in the next Section for the dynamical 
mass model. The position shifts of the orbital centers are rather 
small and can be neglected. The systemic velocity shifts, AVj,,,, 
are computed respect to the nominal value Vjj,j=-300 km s ' . 



Fig. 7. The bottom panel shows the rotation curves obtained for 
the northern and southern halves of the galaxy by averaging the 
data shown in Figure|6]in radial bins 1 kpc wide. Errorbars refer 
to the dispersion around the mean. Filled triangles are for the 
northern half, open squares for the southern half of the galaxy. 
Rotational velocities in the inner region derived from CO data by 
iLoinardet all (119951) are shown with filled and open circles for 
the northern and southern major axis, assumed to be at R4=38° 
and with / = 77°. The filled dots in the upper panel show the 
global rotation curve. As explained in the text, errorbars take 
into account the dispersion around the mean and the differences 
between the values for the two halves of the galaxy. The vertical 
dashed line mark the innermost radius which will be considered 
in this paper for the dynamical analysis. 



3.1. A comparison with otiier rotation curve measures 

In Figure [8] we show for comparison the rotation curves derived 
in this paper (continuous line) and some previously determined 
ones. In panel (a) we display data for north-east galaxy side, in 
panel (b) for the south-west side and in panel (c) the average val- 
ues. Original data in (a) and (b) have been scaled according to 
an assumed distance D = 785 kpc and systemic velocity Vsvs - 
-306 km s ' . We show both optic a l data from.Kent (,1989a^) an d 
21-cm data from lEmers"onl ( Il976h : iNewton & EmersonI (Il977h . 
Unfortunately most of the previous determinations rely on an as- 
sumed fixed inclination for the disk and on data along the major 
axis alone. That implies that possible local velocity perturbations 
will not be averaged out. These are clearly visible especially in 
the literature data between 8 and 20 kpc relative to the south- 
west galaxy half (panel (b)), strongly perturbed by the M32 tidal 
interaction. In the northern side we derive a somewhat lower ro- 
tational velocity, perhaps due to the presence of the warp which 
implies a higher disk inclination, not accounted for by previous 
data analysis. Taking into account these limitations, the general 
agreement seems good. The asteris k symb ols are used in (c) to 
display the average rotation curve of lChem in et al. ( 2009) which 
lies above ours, due especially to somewhat lower inclination 
the authors derive for the tilted ring model. Also, despite their 



analysis masking some perturbations such as the external arm, 
their curve is less smooth than ours. As discussed already in the 
previous Section, this might be due to different choices of ve- 
locity components to extract the rotation curve or to their use of 
weighting function which gives more weight to data points close 
to major axis. A side effect of this choice is to retain any local 
velocity perturbation present along major axis. 

A position velocity plot made along the photometric ma- 
jor axis, at a position angle of 38 °, is shown in Figure |9l Our 
adopted average rotation curve, projected back along major axis, 
has been superimposed to it (diamonds). The average value of 
the disk inclination and systemic velocity, as derived in our best 
fitting tilted ring model in the radial interval of interest, has 
been use d. For a comparison w e display also the average rotation 
curve of lChemin et alj ( l2009h after applying the disk inclination 
and systemic velocity corrections adopted by the authors (aster- 
isk symbols). The figure clearly shows that despite the differ- 
ent gas velocities adopted by the two teams to trace the rotation 
curve, as explained in detail in Section 2, differences in the rota- 
tional velocities at large radii before deprojection are marginal. 
The difference between the two measurements becomes more 
significant when rotation curves (deprojected) are directly com- 
pared. This illustrates the relevance of the differences in the pa- 
rameters of the best tilted ring model found by the two teams. 
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Fig. 8. The rotation curves derived in this paper (continuous 
line) and some previously determined ones. In panel (a) we 
display data for north-east galaxy side, in panel (b) for the 
south-west s ide and in pan el (c) the average values. Filled cir- 
cles are for iKentl (Il989ah optical data along the major axis, 
open squares and open tri angles are for lEmersonl (1 19761) and 
iNewton & Emerso 3 (Il977h 21 -cm data along the major axis. 
Asterisk symbols are used in (c) to display the average rotation 
curve of IChemin et al.l ( l2009l) . Original data in (a) and (b) have 
been scaled according to an assumed distance D = 785 kpc and 
systemic velocity Vsys = -306 km s"'. 



In particular IChemin et al.l (l2009h derive lower inclination an- 
gles for the M3 1 disk and hence higher rotational velocities. The 
agreement between the PV diagram along the major axis and the 
component of the rotational velocities measured along the line 
of sight improves, as it should, when using only data in a very 
small opening angle around the major axis and data in the two 
galaxy halves are kept separate. Notice, for example, the higher 
velocity present along the major axis around 0.7 angular offset 
respect to the average rotation curve. That is also visible in the 
optical data shown in the middle panel of Figure[8] This anoma- 
lous velocity present along the major axis is averaged out when 
data from other directions is taken into account. 

3.2. The bulge-disk decomposition and tlie stellar distribution 

Radial profiles of stellar disk surface brightness can usually 
be represented by an exponential distribution. The mass-to- 
light ratio is determined by the dynamical analysis of the 
rotation curve or by the optical colors once the radial ex- 
ponential disk scale length is known. The disk scale length 
can be established from surface brightness profiles in various 
bands and it might be wavelength dependent. In Andromeda, 
it va ries from 7.7 kpc in the U -band to 5.9 in the R- 
band (IWalterbos & Kennicutill988l) to 4.5 kpc in the K-band 
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Fig. 9. Position velocity diagram along the photometric major 
axis, at a position angle of 38 °. A logarithmic scale is used for 
the brightness and its contours are drawn at 0.013,0.06,0.13,0.6 
Jy/beam. The lowest contour is at the 3-0" level since cr is 
4.2 mJy/beam. The continuous line show the rotation curve 
adopted in this paper w hile the asterisk trace the rotation curve 
of Chemin et al.l ( l2009l) . In both cases, the rotational velocities 
have been corrected using the average value of the disk incli- 
nation and systemic velocity found by the authors in the region 
shown. 

jHiromoto et alj 119831; [Battaner et al.' '198 6|). The im ages used 
at optical wavelengths tWalterbos & Kennicuttlll988i ) can trace 
the disk surface brightness out to galactocentric radii of about 
25 kpc, while the available K-band images (including that ob- 
tained through the 2MASS survey) loose sensitivity beyond 
20 kpc. The discrepancy between the optical and infrared scale- 
length can be easily explained in terms of a radially vary- 
ing star formation history or of an extinction gradient. Light 
in the K-band traces the mass distribution better than optical 
light because of the reduced extinction and because most of 
the stellar mass in galaxy disks is due to old, low mass stars. 
Therefore for the dynamical analysis, K-band scale lengths are 
usually preferred despite the reduced sensitivity of infrared im- 
ages. Recent mid-infrared observations of Andromeda obtained 
with the In frared Array Camera on board of the Spitzer Space 
Telescope (iBarmbv et al.l l2006l) show a disk scale length of 
6.1 kpc at 3.6 //m measured on the same radial interval as mea- 
sured by i Rattaner et al.l (Il986h using the K-band light. A larger 
scale length at 3.6 yum is also found in M33 b y comparing 
the Sp itzer image with the 2MASS K-band image (IVerlev et al.l 
120091) . This might imply that intermediate age, cool supergiants 
contribute a substantial fraction of the NIR emission at 3.6 //m 
(Mould et al. 2008). 

The determination of the disk scale length depends also on 
the bulge-disk light decomposition and K-band images might not 
be deep enough to obtain an accurate bulge-disk decomposition. 
Nevertheless we feel that the shorter disk scalelength measured 
in the K-band images compared to 3.6 fim images is not affected 
by this limitation. M33 in fact has no bulge and indeed the scale 
length at 2.2 /urn is shorter than at 3.6 yum. A spherical model for 
the bulge of M31 is a simplification. Detailed modelling of the 
surface brightness shows that at very least the bulge is an oblate 
spheroid with axis ratio of 0.8 (Kent 1983) but most likely it 
is a triaxial bulge (e.g. lBermanll200ll and references therein). 
For the purpose of this paper, which is the dynamical mass mod- 
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Table 1. The HI rotation curve of M3 1 and the parameters of 
the best fitting tilted ring model 
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elling of the disk rotational velocities (beyond 8 kpc), we can 
neglect the asphericity of the bulge which affects the orbits in 
the innermost few kpc. Using a R'^"^ de Vaucouleurs law (i.e. 
a Sersic index of 4) for the light distribution it is customary 
to characterize the bulge by its effective radius which encloses 
half o f its total light. At optical wav elengths (in the R- and V- 
band) IWalterbos & KennicuttI (Il988h found a bulge effective ra- 
dius of 2.3 kpc and a bulge to disk luminosity ratio B/D=0.82. 
Using the same dataset in the R-band combined with Hubble 
Space Telescope data of Lauer et al. ( 1993) an d data from Kenjj 
(Il983h for the innermost regions IGeehan et al.l (2006) derived a 
much smaller bulge scale radius: 0.61 kpc (similar to what can 
be inferred by inspecting the 2MAS S images). Using the 3.6 jim 
Spitzer image iBarmbv et all ('2006) modelled the bulge light dis- 
tribution with a R'^'* de Vaucouleurs law and found a 1.7 kpc ef- 
fective radius and a bulge to disk light ratio of B/D=0.78. Using 
the same data lSeigar et al.l (l2008l) obtain a bulge effective radius 
of 1.93 kpc for a Sersic index n=1.7, and a disk scalelength of 
5.91 kpc with B/D=0.57. Fits to the bulge light distribu tion us- 
ing sm aller Sersic indexes have also been done bv Worthev et al.l 
(l2005h using a brightness profile from an I-band image out to 
24 kpc: for a Sersic model with n=1.6 the bulge effective ra- 
dius was found to be 0.89 kpc and the disk scale length 5.7 kpc, 
even though a steepening of the scale length is clearly visible be- 
yond 15 kpc. An faint stellar disk which extends as far as 40 kpc 
from the Androm eda nucleus has been recently pointed out by 
llbata et all (l2005h with an exponential scale length of 5.1 kpc in 
the I-band. 

Given the uncertainties in the disk and bulge mass distribu- 
tion we will attempt to fit the rotation curve using a varying disk 



scalelength between 4.5 and 6.1 kpc in steps of 0.2 kpc. For the 
bulge we shall use 4 possible parameter combinations, namely 
an effective radius of 2.0 and 0.7 kpc and a Sersic index n=4, 
and n=1.6. Given the fact that we cannot constrain the dynami- 
cal contribution of the bulge since we are not fitting the motion 
in the inner regions, our purpose will be only to see if the dy- 
namical fit to the disk improves considerably when using any of 
the four combinations. 



3.3. Stellar mass-to-light ratios 

Analysis of the star formation histories of the bulge and disk 
of M31 suggest tha t there is no age difference between the 
bulge and the disk dOlsen et"an l2006l ). However previous at- 
tempts to fit former rotation curves of this galaxy found a higher 
mass-to-light ratio for the disk than for the bulge (Widrow et aT] 
l2003i which was unexpected given the older age of bulges rel- 
ative to dis ks. Hydrodynam ic simulations of the triaxial bulge 
of M3 1 by iBermanI (120011) found a B-band mass-to-light ra- 
tio of 6.5 for the bulge i.e. a stellar mass of 10'*' Mq. In the 
disk of Andromeda there is also a color gradient visible in the 
disk (W alterbos & Kennicutt 1988; Seigar et al. 2008) since B-R 
varies between 1 .7 in the inner regions to 1 .3 in the outer regions 
(this can be due to changes in metallicity, age or extinction). 
According to Bell & de the mass-to-light ratio in the 

K-band expected from B-R colors can vary from the value of 1 in 
the central regions to 0.65 in the outer disk if extinction does not 
change radially. The mass-to-light ratio in the B-band can vary 
between 8 and 2.5 Mq/Lq and we will consider these two values 
as the extreme acceptable disk mass-to-light ratio values. Since 
we fit the rotation curve beyond 8 kpc we cannot constrain the 
:bulge mass-to-light ratio and hence we will consider also models 
with equal mass-to-light ratios for the bulge and the stellar disk. 

The bulge and disk blue luminosities which we shall use 
to compute the mass-to-luminosity ratio are 9 x lO'' and 
2.1 X 10'" Lq respectively. These have b een derived using the in- 
tegrate d B-band magnitude measured bv lWalterbos & KennicuttI 
(1988), corrected for absorption, assuming that the bulge con- 
tribution is 30 % of the total emission in the B-band. As we 
mentioned in the previous subsection, the decomposition of the 
light profile into a bulge and disk component is somewhat uncer- 
tain and especially the bulge integ rated luminosity is not firmly 
estabished yet (see also lKentlll989bl) . 

3.4. The gas surface density 

For the gaseous disk we shall consider the atomic hydrogen and 
the molecular gas surface density. These will be multiplied by 
1.33 to accountfor the presence of helium. Using the best fitting 
tilted ring model (PI) we derive the radial distribution of neutral 
atomic gas, perpendicular to the galactic plane, in the optically 
thin approximation. This is shown in Figure \W\ as a function 
of the galactocentric radius. The corresponding total HI mass 
is 5.4 X 10^ Mq. To consider the possible presence of opaque 
clumps we can multiply the H I surface brightness by 1.3 since 
this is the correction inferred bv Braun et alj (1200 9^. This correc- 
tion has however no noticeable effect for the dynamical analysis 
carried out in the next Section. 

The continuous line in Figure [TO] is the log of the function 
ffjiiR). We shall use fniiR) to compute the dynamical contri- 
bution of the atomic gas mass to the rotation curve. It has the 
following expression 



fmiR) = e 



,2.3(22.0888-0.08759 R) exp(-0.im3 Mp(-0.03275 ij")) 



(6) 
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Fig. 10. The neutral atomic hydrogen column density perpendic- 
ular to the galactic plane is shown as a function of galactocen- 
tric distance. The best fitting tilted ring model PI is used for 
the deconvolution of the 21 -cm line brightness image. The opti- 
cally thin line approximation has been used to convert the sur- 
face brightness in HI gas column density. The continuous line 
shows the total atomic gas surface density of the M31 disk as 
modelled for the dynamical analysis. 

and it provides a good fit to the atomic surface density dis- 
tribution perpendicular to the galactic plane. The sharp decline 
of the HI beyond 27-30 kpc is likely due to the ionization of 
the atomic hydrogen by the local ex tragalactic UV background 
radiation dCorbelli & Salpeterlll993h . Hence the fitting function 
which approximates the total atomic gas distribution is shallower 
than the HI distribution in the outer region. There is also ionized 
atomic gas in t he inner regions, as can be inferred from the Ho- 
emission maps (IWalterbos & BraunI 19941) . However quantifying 
the column density and the geometrical distribution of such ion- 
ized component is subject to several assumptions. The gaseous 
mass involved is small and irrelevant to the rotation curve fit- 
ting since in the inner region it is the stellar component which is 
dynamically dominant. Hence we shall neglect it. 

For the molecular gas fraction we considered the aver- 
age radial variatio n (North-i-South) as shown in Figure 10 of 
lNietenetal.l(l2006l) . 

4. Dynamical analysis 

In this Section we will analyze the mass distribution in the 
Andromeda galaxy. We first attempt to fit the rotation curve of 
Andromeda, derived in the previous Section, from 8 to 37 kpc 
without a dark matter halo. Then we carry on the dynamical anal- 
ysis considering two different models for the radial distribution 
of the dark matter density distribution. We will discuss at the end 
the resulting total mass of Andromeda in the context also of data 
at very large radii from observations of the Andromeda satel- 
lites and other objects. We shall consider both models in which 
the the disk and the bulge have the same mass-to-light ratio and 



models in which these are two independent variables. Unless 
stated differently we shall consider the HI gas in the optically 
thin approximation. The half thickness of the disk is assumed 
to be 0.5 kpc and the contribution of the disk mass component s 
to the rotation curve is computed according to ICasertanol (1 1 9 8 3h . 
We use the reduced chi-square statistic, j^'^, to judge the goodness 
of a model fit. We consider first 2 free parameters: the disk and 
the bulge mass-to-light ratio (M/L)d^b, which we vary continu- 
ously. For the stellar disk we vary the exponential scalelength hd 
using steps of 0.2 kpc in the interval 4.5< < 6.1 kpc. For the 
bulge we consider 2 possible effective radii ht, and Sersic index 
n=4 and n=1.7. 

4.1. Newtonian and non-Newtonian dynamics wittiout dark 
matter 

We first attempt to fit the rotation curve in the framework of 
Newtonian dynamics without considering a dark matter halo. 
The best fit is obtained for a disk scalelength of 6.1 kpc with 
(M/L)d - (M/L)i, - 8.0 Mq/Lq, which is the maximum al- 
lowed value. The fit is slightly better when a bulge effective ra- 
dius of 2 kpc is used. This is shown in the top panel of FigurefTTI 
However the fit is generally poor being the reduced > (!> for 
all possible combinations of parameters (see Table 2). Hence the 
model with no dark matter fails to fit the data under the assump- 
tion of Newtonian gravity. The fit stays poor even if unreason- 
ably high values for the stellar mass-to-light ratio are considered 
(> 8 Mq/Lq). The evident failure is due to the declining rota- 
tion curve predicted by the baryonic mass distribution beyond 
26 kpc. 

An alternative explanation for the mass discrepancy has been 
proposed by Milgrom by means of th e modified Newtonian dy- 
namics or MONO (Milgromllilll. Outside the bulk of the 
mass distribution, MOND predicts a much slower decrease 
of the (effective) gravitational potential, with respect to the 
Newtonian case. This is often suff icient to explain the observe d 
non-keplerian behavior of RCs jSanders & McGaughl l2002h . 
According to this theory the dynamics becomes non-Newtonian 
below a limiting acceleration value, oq ~ 10"** cm s"^, where 
the effective gravitational acceleration takes the value g = 
gnljiix), with gn the acceleration in Newtonian dynamics and 
X - g/oQ. Here we shall use the critical acceleration value 
ao derived from the an alysis of a s ample of rotation curves 
ao = 1.2 X 10"'^ cm s^^ jSanders & M cGaugh 2002). We have 
tested MON D for different choices o f the interpolating function 
fi(x) (see iFamaev & BinnevI I2005L for details). In particular 
we have used the 'standard' and the 'simple' interpolation func- 
tion and found that the 'simple' interpolation function provides 
slightly better fits to the M31 data. Hence we shall use the 'sim- 
ple' interpolating function fj.(x) - x/(l + x) in the rest of the 
paper. 

ICorbelli & Saluccl (l2007h have already tested the former 
M3 1 rotation curve in the MOND framework. They concluded 
that, in M31, MOND fails to fit the falling part of the rotation 
curve at intermediate radii. However, this assessment was made 
using lower quality data and in the absence of an appropriate 
knowledge of the tilted ring model parameters. 

Figure [TT] shows the MOND best-fit model curve when 
{M/L)ij and {M/L)b are two independent free parameters (bot- 
tom panel). The mass-to-light ratio best fitting values are 
(M/L)d = 2.5, (M/L)i, = 5.2, and hd = 4.5 kpc. The value 
of is 1-3, hence MOND provide a good fit to the data. If 
we reduce the number of free parameters to just one by setting 
(M/L)d - (M/L)h the fit is still reasonably good: the lowest is 
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Fig. 11. The 21 -cm rotation curve for M31 and the best fitting 
mass model with no dark halo in the Newtonian case (top panel) 
and for MOND (middle panel). The short-dashed lines indicate 
the Newtonian stellar bulge and disc contribution to the rotation 
curve. The dot-dashed line is the gas contribution. The figure 
displays the the best fitting mass model which is obtained for 
the Newtonian case when (M/L)d - (M/L)b - 8.0 Mq/Lq, - 
6.1 kpc, hh = 2.0 kpc and n=1.6. The MOND best fit (x^ = 1.3) 
is obtained when the free parameters are = 4.5 kpc, (M/L)d - 
2.5 Mq/Lq (M/L)t = 5.2 Mq/Lq, and hi, = 2 kpc, n=1.6. For 
comparison we show in the bottom panel the MOND best fit 
when the disk scalelength is 6.1 kpc. The reduced in this 
case is higher (x^ = 3.1), (M/L)d = 4.1 Mq/Lq, and {M/L)b = 
3.1 Mo/Lo . 



1.37 coiTesponding to (M/L)d = (M/L)t = 3.3 Mq/Lq. The ro- 
tational velocities predicted by MOND are only slightly higher 
than observed for 20 < /? < 30 kpc and slightly lower than the 
data for 10 < < 20 kpc. As shown in Table 2, the goodness 
of the fits are not very sensible to the bulge mass-to-light ratio, 
since we are excluding the central regions. It is however to be 
noticed that MOND best fits to the actual data require a higher 
mass-to-light ratio for the bulge than for the disk, in agreement 
with the expected older age of the bulge stellar component. In 
M31 it is around 10 kpc that g„ ~ ao and non-Newtonian cor- 
rections start to be important and force a falling Newtonian RC 
into a flat one, more consistent with the data. We notice that for 
MOND the value of the disk scalelength used to fit the bary- 
onic matter distribution is very important. In fact the fit becomes 
poor if /zj/ = 6.1 kpc, as shown in Figure [TT] If future photomet- 
ric studies will confirm a disk scalelength of 6.1 kpc then one 
will have to consider possible variations of the assumed distance 
to M31 to make MOND predictions more consistent with the 
kinematics traced by 21 -cm data. 

4.2. Dark matter halo models 

In the previous subsection we have seen that Newtonian dynamic 
fits to the rotation curve without considering a dark matter halo 
are rather poor We will now use the M31 rotation curve pre- 
sented in this paper to test in detail the consistency of a possi- 
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Fig. 12. The M3 1 rotation curve (points) and the best-fitting 
mass models (solid line) using a Burkert dark halo profile with 
hd = 5.1 kpc, hh = 2 kpc and n=4. Also shown are the dark halo 
contribution (dot-dashed line), the stellar disk and bulge (short- 
dashed line) and the gas contribution (long-dashed line). In the 
top panel, we show the best fit mass model (x^ - 0.81) with 
(M/L)b = 4.5 Mo/Lo, (M/L)d = 8.0 MJLq and Rb = 77 kpc. 
The case shown in the bottom panel refers to a fixed, lower value 
of the core radius, namely Rb = 28 kpc. For this case the best 
fitting values of the mass-to-light ratios are {M/L)h = 4.9 Mq/L© 
and (M/L)d = 7.4 MJL^ (x^ = 1.17). 



ble halo density profile with theoretical models which predict a 
well defined dark matter distribution around galaxies. Namely, 
we shall consider a spherical halo with a da r k mat t er den sity 
profile as originally derived by iNavarro et al. 1 (119961 Il997h for 
galaxies forming in a Cold Dark Matter scena rio. We co nsider 
also the Burkert dark matter density profile (Burkertl ll995^ since 
this successfully fitted the rotation curve of da rk matter domi- 
nated dwarf galaxies (e.g. i Gentile et al.l 120071 and references 
therein). Both models can describe the dark matter halo density 
profile u sing two parameters. The density profile proposed by 
iBurkertI (fl995i) has a constant-density core and is given by: 



p{R) 



Pb 



(7) 



(1 + 



A strong correlation between the two parameters pB and Rb has 
been found by fitting the r otation curve of low mass disk galaxies 
dSalucci & BurkertlboOOl) . namely 



. Rb ^'^ 
Mb =4.3x10^—^ Mo 
kpc 



with Mb = \ .6pbR\ 



(8) 



Using this relation, we show in Figure [T2] the best fitting 
mass model when a Burkert model for the dark halo is con- 
sidered. The values of the free parameters are: Rb - 77 kpc, 
(M/L)d - 8.0 Mq/Lq and hd - 5.1 kpc. The reduced value 
is close to unity, = 0.81, meaning the fit is generally good. 
The best fit does not depend much on the assumed bulge pa- 
rameters. Even considering {M/L)d - {M/L)b the halo model 
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Fig. 13. The 68 % (darker regions) and 95.4 % (lighter regions) 
confidence areas in the {M/Q^-Rg plane. Clearly the confidence 
areas extend even beyond the largest Rg value shown. In panel 
(a) (M/L)d = iM/L)h, while in (b) {M/L)d and (M/Lh are two 
independent variables. 

provides an excellent fit to the data. We define the virial mass, 
Mg'', in the case of a Burkert halo as the mass enclosed in the 
sphere which has the average dark matter density equal to 98 pcr 
where the critical density for closure is p„ = 3H^/8nG with 
Hq = 71 km s The virial mass of the dark halo for best fit- 
ting mass model is uncomfortably high: M^" = 1.2 10'^ Mq. A 
Burkert halo which fits M3 1 is indistinguishable from a constant 
density dark matter model since even at the outermost fitted ra- 
dius, the dark matter density does not decline yet as an R^^ or 
R^^ power law. If we consider core radii smaller than the best 
fitting value, comparable to the extent of the HI disk, Burkert 
halo models still provides acceptable fits to the data and implies 
lower virial masses. In Figure[T3]we show the 68 % and 95.4 % 
confidence areas in the {M/L)ci-Rb plane. For 22 < Rb < 37 kpc 
we have virial masses 8.3 lO" < M'j^' < 2.5 lO'^ Mq and;^'^ 
values in the 95.4 % confidence area (x^ increases as Rg de- 
creases). In the bottom panel of Figure[T2]we show, as an exam- 
ple, a similar mass model to that shown in the top panel, except 
for the value of the core radius which we set equal to 28 kpc. 
With such low value of Rb the Burkert halo model has a virial 
mass of M™' =1.4 10^^ M© and provides still an acceptable fit 
to the data Or^ = 1.17). 

The NFW density profile is: 



P(R) 



Pnfw 



Rnfw \ Rnfw I 



(9) 



Numerical simulations of galaxy formation find a correlation 
between pnfw a nd Rnfw which depends on the cosmologi- 
cal model (e.g . iNavarro et all 119971: lAvila-Reeseet all 120011: 
lEke et al.l200UlBullock et al.l200 1'). Often this correlation is ex- 
pressed using the concentration parameter C = Ra/Rnfw and 
Ma or Va. Ra is the radius of a sphere containing a mean den- 
sity A times the cosmological critical density, Va and Ma are the 



characteristic velocity and mass inside R^. In this paper we use 
the results of N-body simulations in ACDM, cosmology to de- 
fine the relation between the concentration and the virial mass 
of dark matter halos. We adopt a fla t ACDM cosmolog y with 
parameters from the WMAP results (iSpergel et al.ll2003l) . mat- 
ter density Q.m = 0.27, baryon density Q./, - 0.044, a normal- 
ized Hubble constant h = //()/(100 km s"' Mpc"')=0.71. In this 
framework, we use the results of the numerical simulation of 
Maccio et al. (2007), which are in agreement with the results of 
Neto et al., (2007,) . and give the following relation between the 
concentration and the halo mass for relaxed halos: 



C 



7.5(Ma/101VMo)-°°''« 



(10) 



where Ma is the halo mass for a density contrast A = 98. 
Generally the dispersion in the concentration decreases mono- 
tonically as a function of halo mass. The scatter around logC 
given by the above relation, is about +0.13 at the expected virial 
mass of Andromeda (i.e. 10^^ M©), and a similar value is found 
for the average total scatter over the mass range considered. 

We now fit the M3 1 rotation curve using two free parame- 
ters: the halo concentration C and (M/L)^, the disk mass-to-light 
ratio. We shall consider first no dispersion around the relation 
given by the above equation and a bulge mass-to-light ratio equal 
to that of the disk. The 1-cr, 2-cr and 3-cr ranges are determined 
using the reduced and assuming Gaussian statistics. They are 
computed by determining in the free parameter space the projec- 
tion ranges, along each axis, of the hypersurfaces corresponding 
to the 68.3% and 95.4% confidence levels. 

The best fitting mass model is obtained for a disk scale- 
length hci = 4.5 kpc. The bulge parameters are not of much 
relevance for the goodness of the fit in the region of interest 
to this paper The stellar mass-to-light ratio for the best fit is 
(M/L)d - 4.2 Mq/Lq, if we assume that it does not vary between 
the bulge and the disk component, and the value of the reduced 

is 1.12. The total dark halo mass is Ma = 10'^ M© (corre- 
sponding to C = 11.9). If we allow variations between the disk 
and the bulge mass-to-light ratio, the best fit mass model gives 
(M/L)d - 5.0 Mq/Lq and (M/L)b - 2.1 Mq/Lq and a minimum 

value of 1.08. These combination of M/L ratios is not real- 
istic since we don't expect the bulge to have a lower M/L ratio 
than the disk. But higher values of the bulge M/L ratio increases 
only slightly the x^ value. In Figure [14] we show the modelled 
rotation curve according to the best fitting mass model when 
hb - 2.0 kpc and n=4. The dai'k halo mass is Ma = 1 .2x10'^ Mq 
corresponding to C = 12. As we increase the disk scalelength 
the x^ increases slightly, the dark halo mass decreases and the 
stellar mass-to-light ratio increases (for example for a change in 
the disk scalelength from 4.5 to 6.1 kpc the minimum;t'^ changes 
from 1.12 to 1.34, Ma decreases from 1.2 x lO'^ to 7.5 x lO" Mq 
and {MIL)h4 = 5.7 M^/Lo). 

We have computed the confidence area in the (M/L)ij - C 
plane and shown them in Figure[T5] When the mass-to-light ratio 
of the disk and the bulge are two independent variables the best 
fitting (M/L)i value is unrealistically low but the value of the 
concentration parameters does not change. FigurefTSlshows also 
the wider confidence areas obtained when twice the dispersion 
of ±0.13, as estimated from numerical simulations, is considered 
around logC in the C-Ma relation. 

Table |2] summarizes the main results of the mass models 
considered for fitting the rotation curve of M31. In column (1) 
we give the short name of the mass model: No-DM is simply 
a Newtonian dynamic fit with no dark matter For MOND we 
use the no dark matter and the "simple" interpolation function. 
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Table 2. Parameters of the best fitting mass model. 



MODEL 




hi, 




(M/L),_, 


(M/L)rf 




T 






kpc 


kpc 


(M/L)o 


(M/L)o 






No-DM 


4 


0.7 


6.1 


8.0 


8.0 




6.10 


No-DM 


4 


2.0 


6.1 


8.0 


8.0 




6.15 


No-DM 


1.6 


0.7 


6.1 


8.0 


8.0 




6.17 


No-DM 


1.6 


2.0 


6.1 


8.0 


8.0 




6.02 


MOND 


4 


0.7 


4.5 


5.2 


2.5 




1.32 


MOND 


4 


2.0 


4.5 


5.6 


2.5 




1.41 


MOND 


1.6 


0.7 


4.5 


5.1 


2.5 




1.31 


MOND 


1.6 


2.0 


4.5 


5.2 


2.5 




1.30 


Burkert 


4 


0.7 


5.3 


4.6 


8.0 


Rb=S3 kpc 


0.84 


Burkert 


4 


2.0 


5.1 


4.5 


8.0 


i?s=77 kpc 


0.81 


Burkert 


1.6 


0.7 


5.3 


4.5 


8.0 


Rb=74 kpc 


0.84 


Burkert 


1.6 


2.0 


5.3 


4.5 


8.0 


Rb=S3 kpc 


0.83 


ACDM 


4 


0.7 


4.5 


2.6 


4.9 


C=11.9 


1.09 


ACDM 


4 


2.0 


4.5 


2.7 


5.0 


C=12.0 


1.08 


ACDM 


1.6 


0.7 


4.5 


2.6 


4.9 


C=11.9 


1.10 


ACDM 


1.6 


2.0 


4.5 


2.6 


4.9 


C=11.9 


1.08 


MOND 


4 


0.7 


4.5 


3.3 


3.3 




1.37 


MOND 


4 


2.0 


4.5 


3.3 


3.3 




1.48 


MOND 


1.6 


0.7 


4.5 


3.3 


3.3 




1.37 


MOND 


1.6 


2.0 


4.5 


3.3 


3.3 




1.37 


Burkert 


4 


0.7 


5.5 


7.0 


7.0 


Rb=74 kpc 


0.90 


Burkert 


4 


2.0 


5.3 


7.0 


7.0 


i?fl=78 kpc 


0.85 


Burkert 


1.6 


0.7 


5.5 


6.9 


6.9 


Rb=66 kpc 


0.92 


Burkert 


1.6 


2.0 


5.5 


7.0 


7.0 


Rb=16 kpc 


0.89 


ACDM 


4 


0.7 


4.5 


4.1 


4.1 


C=11.8 


1.14 


ACDM 


4 


2.0 


4.5 


4.2 


4.2 


C=11.9 


1.12 


ACDM 


1.6 


0.7 


4.5 


4.0 


4.0 


C=11.8 


1.15 


ACDM 


1.6 


2.0 


4.5 


4.1 


4.1 


C=11.8 


1.12 



The scalelengths of the bulge and of the disc are labeled with the 
symbol hh and respectively. We consider h^ as a free parame- 
ter in the interval 4.5-6. 1 kpc. The best fitting value of a possible 
additional free parameter of each model, P+ is given in column 
(7). For the No-DM and MOND models there is no additional 
free parameter. For the Burkert halo model P+ is the core radius 
in kpc, for the ACDM mass model P^ is the concentration pa- 
rameter C. For the bulge and disk mass-to-light ratio, the range 
of possible values considered are: 2.5<(M/L)/, < 8.0 in units of 
(M/L)0. In the first half of the Table, we consider both (M/L)i 
and (M/L)rf as free parameters while in the second part we use 
(MfL)b=(M/L)d. Since, for the simple Newtonian dynamic fit 
with no dark matter, the bulge and disk mass-to-light ratios come 
out equal to the highest allowed value, the parameters are listed 
only in the first part of the Table. 

We can see clearly in Table|2]that the goodness of the fit does 
not depend on the bulge mass distribution, as expected, since we 
are only fitting the rotation curve for > 8 kpc. It also shows that 
only a pure Newtonian mass model with no dark matter halo fails 
to fit the data. Both MOND, a Burkert halo model with a large 
constant density core, and the NFW profile in the framework of 
ACDM provide good mass model fits to the rotation curve of 
Andromeda derived in this paper In the next subsection we shall 
use data from other sources to test the rotation curve at larger 
radii than those provided by the HI dataset. 

4.3. Baryonic and total mass of Andromeda 

We consider now the mass estimate of Andromeda from differ- 
ent sources at galactocentric radii from 30 to 560 kpc. We de- 
rive a rotation curve by computing the expected circular rota- 
tional velocity given the mass estimate within a radius R. We 



consider data from planetary nebulae, globular clusters, stellar 
streams and Andromeda satellites to constrain the Andromeda 
total mass at large radii. Since each paper considers an ensemble 
of objects and describes in detail the method used to derive the 
mass, we give in Table 3 the resulting mass estimate and refer- 
ence the original papers where a more detailed description of the 
database and analysis can be found. If the authors have assumed 
a different distance to Andromeda than what we use in this pa- 
per, we have scaled their radius and mass estimate accordingly. 
In Figure [16] we plot the rotational velocities at large distances 
from the center of Andromeda derived from the mass estimates 
given in Table 3. 

We plot also the velocities at large radii predicted by the 
ACDM mass models which best fits our rotation curve i.e. a 
NFW dark halo profile with concentration C = 12 (continu- 
ous line). A Burkert halo model with a very large core radius, 
such as that required by the best fitting mass model to the 21- 
cm rotation curve data, predicts higher velocities than observed 
at distances shown in Figure [16] But a Burkert halo profile with 
a core radius iJ^ = 28 kpc, compatible with the 21 -cm rotation 
curve data, is also compatible with data at large galactocentric 
distances (short dashed line in Figure [T6]l. For both halo models 
we have considered the values of the free parameters hd, (M/L)d 
and {M/L)t,, which best fit 21 -cm rotation curve. Given the fact 
that the ACDM and Burkert halo were constrained using only 
HI data between 8 and 37 kpc it is quite remarkable how well 
the two halo models fit the data out to about 560 kpc. We did not 
test MOND predictions in this framework because Andromeda 
masses in the original paper have been derived using Newtonian 
dynamic. 

Since the best fit Burkert halo model was essentially a con- 
stant density profile in the region traced by the 21 -cm rotation 
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Fig. 14. The M3 1 rotation curve (points) and the best-fitting 
mass model (solid line) using the NFW dark halo profile with 
C - 12 in the frame of AC DM. Also shown are the dark halo 
contribution (dot-dashed line), the stellar disk and bulge (short- 
dashed line) and the gas contribution (long-dashed line). The 
bottom panel refers to the case ((M/L)^ = iM/L)i, = 4.2M0/L0). 
The top panel refers to the best fit when the mass-to-light ratio 
of the disk and the bulge are two independent variables. For the 
best fit {M/L)d = 5.0 and an {M/L)b = 2.7 Mq/Lq. These com- 
bination of M/L ratios is not realistic since we don't expect the 
bulge to have a lower M/L ratio than the disk. But higher values 
of the bulge M/L ratio increases only slightly the value and 
require similar dark matter distribution. Both fits refer to a bulge 
model with hh = 2.0 kpc and n=4. 

curve, we show also the predicted velocities of this model in 
case the dark halo is truncated beyond 37 kpc. A simple keple- 
rian fall off of the observed velocities outside the region covered 
by our 21 -cm dataset however fails to fit the data at larger galac- 
tocentric radii. Only the outermost 3 data points are compatible 
with a keplerian fall-off regime. For the KCDM mass model, the 
virial radius corresponding to a concentration C - 12 is 270 kpc, 
hence the predicted rotational velocities of the outermost data 
points in Figure [16] are in the keplerian regime. 

The best fitting ACDM mass model implies a stellar 
(disk-Hbulge) mass of 1 .3 x 10" Mq and a dark matter halo mass 
of 1.2 X 10'^ Mq. Adding to the stellar mass the cold gas mass 
(neutral and molecular hydrogen plus helium) of 7.7x10'' M© we 
estimate a total baryonic mass of Andromeda of 1.4 x 10" Mq. 
This adds to the dark matter halo mass, giving 

M,o, = 1.3x lO'^Mo (11) 

as our best estimate of the total mass of the Andromeda 
galaxy. The associated baryonic fraction is 0.12, very similar to 
the cosmic inferred value of 0.14. 

5. Summary 

We determine the rotation curve of the M31 disk from 8 to 
37 kpc using a tilted ring model fit to the velocity field mapped 



Fig. 15. The 68 % (darker regions) and 95.4 % (lighter regions) 
confidence areas in the {M/L)^-C plane. Panel (a) and (b) refers 
to the case where the average C-Ma relation is considered with 
no dispersions around it. In panel (a) (M/L)d - {M/L)b, while in 

(b) (M/L)d and {M/L)h are two independent variables. In panel 

(c) the confidence areas are for the {M/L)d - (M/L)h case when 
2-0" = 0.26 around the average logC-logMA relation is consid- 
ered. 

in the full-disk, 21 -cm imaging survey of lBraun et all ( l2009l) . 
The orientation of the rings have been determined using three 
different techniques which give rather similar results. For our 
most accurate modelling method (PI), we use 1 1 equally spaced 
free rings, which cover galactocentric distances between 9 and 
36 kpc, whose parameters are varied independently. Each free 
ring has 6 degrees of freedom, since we allowed the systemic 
velocity and center position of each ring to vary (in addition to 
the circular rotation, inclination and position angle). This im- 
plies a total of 66 degrees of freedom in our model. Between 
two consecutive free rings, parameters are determined by a linear 
interpolation. We find that the disk of M3 1 warps from 25 kpc 
outwards by lowering its position angle and becoming more in- 
clined with respect to our line of sight. The disk reaches an in- 
clination of 86° at 35 kpc. The geometry of the outermost two 
rings has somewhat larger uncertainties, but the tilted ring model 
which gives the best fit to the data also produces consistent rota- 
tion curves in the two separate halves of the galaxy. Furthermore, 
these rotation curves do not depend on the value of limiting an- 
gle around the major axis chosen for selecting the data. We find 
-306 km s"' as the average value of the systemic velocity of 
the gaseous disk of M3 1 . The rotation curve of M3 1 is consis- 
tent with being flat beyond 20 kpc and we carry on a dynamical 
analysis to determine the baryonic and non-baryonic mass dis- 
tribution of the nearest spiral galaxy. 

The M31 rotation curve cannot be reproduced using 
Newtonian dynamics and only the stellar and gaseous mass com- 
ponents. Without a dark matter halo however, MOND provide a 
good fit to the galaxy gravitational potential in the region consid- 
ered. We test the density profile and mass predictions of hierar- 
chical clustering and structure formation in a ACDM cosmology, 
together with a dark halo model having a constant density core. 
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Table 3. Andromeda mass estimated at large galactocentric radii. 



Reference 


R 


Mass 




Objects 




[kpc] 


[10"' Mo] 


1 -1 

km s 




Evans & Wilkinson (2000) 


32 


■-'^-10 


230+= 


17 Globular Clusters,^ 


This paper 


37 


49+12 


240^1 


21 -cm data 


Evans & Wilkinson (2000) 


41 


4C+34 


99^+69 

-63 
9ne+/ 
208_g 


17 Globular Clusters,^ 


T f^t q1 /^onns^ 

Lcc ct ai. (ZUUoj 


CC 

JJ 


cc+4 
J J ^ 


jU'4- OiODUloT V^lUSlciS 


Galleti et al. (2006) 


60 


44+26 


178+f 


349 Globular Clusters 


Cote et al. (2000) 


100 


79+5 
75+25 

-13 


185+« 


12 Satellites,^ 


Ibata et al. (2004) 


125 


161+25 


Stellar Stream 


Fardal et al. (2006) 


125 


74+12 


160!jl 


Stellar Stream 


Courteau & van den Bergh (1999) 


268 


137+|« 


149+9 
98!« 


7 Satellites 


_Evans & Wilkinson (2000) 


560 


1 9C+180 
^^-^-60 


11 Satellites^ 


Evans et al. (200Q) 


560 


QQ+146 

^^-63 


87!5' 


16 Satellites^ 



" Best mass estimate using GCs; the lower mass extreme is the value inferred using 9 PNe, the upper mass extreme is from data on 1 1 Satellites 

* Uncertanties derived from "Error analysis" section of original paper 
Uncertanties include uncertanties on the isotropy degree of the velocity distribution and the mass estimate using only the 14 most distant GCs 
Mass estimate assuming isotropic orbits; the mass range is larger, 19*_™, if models with only radial or circular orbits are included 
Uncertanties derived from "Error analysis" section of original paper 

f Mass estimate assuming Osipkov-Merrit distribution function; the given mass range includes also the mass estimated assuming a constant 
orbital anisotropy. 
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Fig. 16. Rotational velocities predicted at large distances from 
the center of Andromeda according to several datasets analyzed 
by previous papers (see Table 3 for details). The continuous line 
shows the rotational velocities predicted by an extrapolation at 
large radii of the best fitting mass model with a NFW dark halo 
profile (C - 12). The short dashed line shows the predicted ve- 
locities of a Burkert halo model with Rb - 28 kpc. The long 
dashed line is for a constant halo dark matter density profile 
which is truncated outside 37 kpc and gives a Keplerian fall off 
of the velocity at larger distances. 



fits M3 1 has a core radius comparable to the size of the disk of 
M3 1 and therefore is in practice a constant dark matter density 
model. 

Using the relation between the concentration parameter C 
and the dark halo mass Ma as for a NFW density profile in 
a flat ACDM cosmology, we find a best fit halo concentra- 
tion parameter C=12 which implies a dark matter halo mass 
Ma= 1.2x10'^ Mq. If we assume that the stellar disk and the 
bulge have the same mass-to-light ratio we find (M/L)ij = 
(M/L)b - 4.2 ± 0.4 Mq/Lq. If the mass-to-light ratio of the 
disk and the bulge are two independent variables then the best 
fit gives a slightly higher value for the disk, (M/L)d - 5.0+j q. 
We are unable to constrain the bulge mass-to-light ratio since we 
discarded the innermost rotation curve in our fit (derived with- 
out considering elliptical orbits). A wider range of C and {M/L)d 
values are found when a dispersion is considered around the av- 
erage log C - log Ma relation, as suggested by numerical simu- 
lations of structure formation in a ACDM cosmology. 

An interesting result is that the extrapolation of a constant 
core dark halo model, as well as of the best fit ACDM dark halo 
model, beyond the region traced by the rotation curve are in good 
agreement with the Andromeda mass traced by other dynamical 
indicators (globular clusters, streams, satellites) at larger radii, 
out to 560 kpc. The constant-core best fitting halo model has a 
very large core radius (77 kpc) and a high virial mass, not con- 
sistent with the data at large galactocentric radii. However mod- 
els with a somewhat smaller core radius, provide acceptable fit 
to the 21 -cm rotation curve and to data at larger galactocentric 
radii. The total estimated mass of M3 1 from our mass model fit 
to the 21 -cm rotation curve in the framework of ACDM cos- 
mology is 1.3+y3 X 10'^ Mq, with a 0.12 baryonic fraction. This 
is similar to the cosmic inferred baryonic fraction of 0.14 and 
implies a formation redshift Zf - 1.2 for the Andromeda galaxy. 
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Both models are able to reproduce the rotation curve of M3 1 to 
a high level of accuracy. The constant density core model which 
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